#library(ncdf)

##Funktion mit den Laengen-und Breitengraden als Eingabeparameter
##x ist die Eingabematrix
cut_matrix <- function(x,
                       Res,   
                       LongMin, 
                       LongMax, 
                       LatMin, 
                       LatMax,
                       cutLongMin,
                       cutLongMax,
                       cutLatMin,
                       cutLatMax) {
  
  ## calculate row and column names
  ## need to be in the pixel centre (i.e. +/- 0.5 * Res) 
  ## otherwise we get one too many entries in the seq() command
  colseqmin <- LongMin + 0.5 * Res
  colseqmax <- LongMax - 0.5 * Res
  rowseqmin <- LatMin + 0.5 * Res
  rowseqmax <- LatMax - 0.5 * Res
  ##mit seq() Datensequenz erstellen von colseqmin bis colseqmax mit den Abst??nden Res
  colnms <- seq(colseqmin, colseqmax, Res)
  rownms <- seq(rowseqmin, rowseqmax, Res)
  
  ## set row and col names so that they mirror what is supplied by LongMin etc.
  ##mit dimnames() lassen sich den Dimensionen einer Matrix "Labels" hinterlegen
  dimnames(x) <- list(rownms, colnms)
  
  ## same as above to calculate region to be cut
  cutlatmin <- cutLatMin + 0.5 * Res
  cutlatmax <- cutLatMax - 0.5 * Res
  cutlonmin <- cutLongMin + 0.5 * Res
  cutlonmax <- cutLongMax - 0.5 * Res
  
  ## now we extract which() indeces in row and col names correspond to 
  ## supplied cut points for lat and lon
  fromlat <- which(dimnames(x)[[1]] == cutlatmin)
  tolat <- which(dimnames(x)[[1]] == cutlatmax)
  fromlon <- which(dimnames(x)[[2]] == cutlonmin)
  tolon <- which(dimnames(x)[[2]] == cutlonmax)
  
  ## fnally, we extract respective rows and columns identified above
  x_matrix <- x[fromlat:tolat, fromlon:tolon]
  
  x_matrix
  
}
  
  
  
#   ##Berechnung der Laengengeraden  
#   LatitudeMin <- ifelse(LatMin < 0, (LatMin+180)/Res+1, LatMin/Res+720)
#   LatitudeMax <- ifelse(LatMax < 0, (LatMax+180)/Res+1, LatMax/Res+720)
#   
#   ##Berechnung der Breitengraden
#   LongitudeMin <- ifelse(LongMin < 0, (LongMin+50)/Res+1, LongMin/0.25+200)
#   LongitudeMax <- ifelse(LongMax < 0, (LongMax+50)/Res+1, LongMax/0.25+200)
#   
#   ##wenn man den combine Befehl ausfuehrt, gibt er die gesuchten Spalten- und 
#   ##Zeilengrenzen aus
#   ##prec_cut <- c(LatitudeMin, LatitudeMax, LongitudeMin, LongitudeMax)
#   
#   ##aber das direkt zurechtschneiden funktioniert noch nicht =(
#   x_matrix <- x[(LatitudeMin:LatitudeMax), (LongitudeMin:LongitudeMax)]
#   
#   x_matrix
# 
#   }

##Daten einladen und Niederschlagsvariable (= Eingabematrix) extrahieren
#data <- open.ncdf("C:/Users/Steffi/Documents/=Uni/Projektarbeit/BspProjektarbeit.nc", write = FALSE, readumlim = TRUE)
#prec <- get.var.ncdf(data, varid = "precipitation")

##Funktion ausprobieren
##Bei der Eingabe der Parameter muss beachtet werden, dass LongMin < LongMax und 
##LatMin < LatMax ist, ausser wenn der Ausschnitt nach Osten ueber den -180sten 
##Laengengrad hinausgeht


## EXAMPLE 1 ##

## define vector of 4 * 180 * 360 random values (corrsponding to 0.25 deg. res)
v <- rnorm(180*360*4, 70, 40)
## write values of v into matrix
x <- matrix(v, nrow = 180*4, ncol = 360*4)

## cut marix with function
Ausgabe <- cut_matrix(x = x,
                      Res = 0.25,
                      LongMin = -180,
                      LongMax = 180,
                      LatMin = -90,
                      LatMax = 90,
                      cutLongMin = 36,
                      cutLongMax = 38,
                      cutLatMin = -4,
                      cutLatMax = -2)

levelplot(t(Ausgabe))

## EXAMPLE 2 ##
## define vector of 100 * 360 random values (corresponding to TRMM at 1 deg. res)
w <- rnorm(180*360, 70, 40)
## write values of v into matrix
y <- matrix(w, nrow = 100, ncol = 360)

## cut marix with function (this time reference coordinate system from 0 to 360)
Ausgabe2 <- cut_matrix(x = y,
                       Res = 1,
                       LongMin = 0,
                       LongMax = 360,
                       LatMin = -50,
                       LatMax = 50,
                       cutLongMin = 0,
                       cutLongMax = 20,
                       cutLatMin = 13,
                       cutLatMax = 28)

levelplot(t(Ausgabe2), col.regions = colorRampPalette(c("white", "blue", "red")))
